GEOG-533 Lab 9

Question 1

Use R package UScensus2010county to complete the following tasks: (20 pt.)

Question 1(a)

Plot a map of New York counties using the plot function.

if(!require('UScensus2010')) install.packages('UScensus2010')
library(UScensus2010)
if(!require('UScensus2010county')) install.county('osx')
library(UScensus2010county)
data('new_york.county10')
shp <- new_york.county10
plot(shp)

Question 1(b)

Plot a map of New York counties using the qtm function.

if(!require('tmap')) install.packages('tmap')
library(tmap)
qtm(shp,fill = 'P0010001')

Question 1(c)

How many counties in New York State?

df <- shp@data
df

nrow(df)

[1] 62

Question 1(d)

What’s the 3-digit fips code of Broome County?

str(df)

'data.frame':   62 obs. of  462 variables:
 $ NAME10    : chr  "Broome" "Seneca" "Schenectady" "Niagara" ...
 $ NAMELSAD10: chr  "Broome County" "Seneca County" "Schenectady County" "Niagara County" ...
 $ state     : chr  "36" "36" "36" "36" ...
 $ county    : chr  "007" "099" "093" "063" ...
 $ fips      : chr  "36007" "36099" "36093" "36063" ...
 $ P0010001  : num  200600 35251 154727 216469 159429 ...
 $ P0020001  : num  200600 35251 154727 216469 159429 ...
 $ P0020002  : num  0 0 0 0 0 0 0 0 0 0 ...
 $ P0020003  : num  0 0 0 0 0 0 0 0 0 0 ...
 $ P0020004  : num  0 0 0 0 0 0 0 0 0 0 ...
 $ P0020005  : num  0 0 0 0 0 0 0 0 0 0 ...
 $ P0020006  : num  200600 35251 154727 216469 159429 ...
 $ P0030001  : num  200600 35251 154727 216469 159429 ...
 $ P0030002  : num  176444 32591 123211 191673 139529 ...
 $ P0030003  : num  9614 1607 14710 14851 10338 ...
 $ P0030004  : num  396 104 575 2285 385 ...
 $ P0030005  : num  7065 244 4960 1823 3517 ...
 $ P0030006  : num  82 2 105 62 34 34 213 16 8 19 ...
 $ P0030007  : num  1912 240 5364 1108 1783 ...
 $ P0030008  : num  5087 463 5802 4667 3843 ...
 $ P0040001  : num  200600 35251 154727 216469 159429 ...
 $ P0040002  : num  193822 34299 145900 211775 153349 ...
 $ P0040003  : num  6778 952 8827 4694 6080 ...
 $ P0050001  : num  200600 35251 154727 216469 159429 ...
 $ P0050002  : num  193822 34299 145900 211775 153349 ...
 $ P0050003  : num  173074 31999 119409 188907 136555 ...
 $ P0050004  : num  8850 1513 13528 14511 9592 ...
 $ P0050005  : num  328 96 445 2135 283 ...
 $ P0050006  : num  7019 238 4917 1807 3469 ...
 $ P0050007  : num  60 2 81 55 25 29 137 14 5 15 ...
 $ P0050008  : num  242 47 2708 185 241 ...
 $ P0050009  : num  4249 404 4812 4175 3184 ...
 $ P0050010  : num  6778 952 8827 4694 6080 ...
 $ P0050011  : num  3370 592 3802 2766 2974 ...
 $ P0050012  : num  764 94 1182 340 746 ...
 $ P0050013  : num  68 8 130 150 102 183 1000 55 27 20 ...
 $ P0050014  : num  46 6 43 16 48 46 467 14 7 8 ...
 $ P0050015  : num  22 0 24 7 9 5 76 2 3 4 ...
 $ P0050016  : num  1670 193 2656 923 1542 ...
 $ P0050017  : num  838 59 990 492 659 ...
 $ P0060001  : num  206107 35732 161052 221479 163567 ...
 $ P0060002  : num  181009 33034 127516 195911 143012 ...
 $ P0060003  : num  12499 1812 17992 17632 12620 ...
 $ P0060004  : num  1664 281 1769 3859 1273 ...
 $ P0060005  : num  8193 291 6587 2407 4293 ...
 $ P0060006  : num  243 8 391 158 123 165 683 65 34 58 ...
 $ P0060007  : num  2499 306 6797 1512 2246 ...
 $ P0070001  : num  206107 35732 161052 221479 163567 ...
 $ P0070002  : num  198382 34718 151125 216232 156746 ...
 $ P0070003  : num  176969 32386 123004 192737 139499 ...
 $ P0070004  : num  11278 1705 16273 17058 11487 ...
 $ P0070005  : num  1445 269 1465 3589 1058 ...
 $ P0070006  : num  8072 285 6436 2354 4205 ...
 $ P0070007  : num  191 8 327 141 100 122 439 51 25 42 ...
 $ P0070008  : num  427 65 3620 353 397 ...
 $ P0070009  : num  7725 1014 9927 5247 6821 ...
 $ P0070010  : num  4040 648 4512 3174 3513 ...
 $ P0070011  : num  1221 107 1719 574 1133 ...
 $ P0070012  : num  219 12 304 270 215 ...
 $ P0070013  : num  121 6 151 53 88 ...
 $ P0070014  : num  52 0 64 17 23 43 244 14 9 16 ...
 $ P0070015  : num  2072 241 3177 1159 1849 ...
 $ P0080001  : num  200600 35251 154727 216469 159429 ...
 $ P0080002  : num  195513 34788 148925 211802 155586 ...
 $ P0080003  : num  176444 32591 123211 191673 139529 ...
 $ P0080004  : num  9614 1607 14710 14851 10338 ...
 $ P0080005  : num  396 104 575 2285 385 ...
 $ P0080006  : num  7065 244 4960 1823 3517 ...
 $ P0080007  : num  82 2 105 62 34 34 213 16 8 19 ...
 $ P0080008  : num  1912 240 5364 1108 1783 ...
 $ P0080009  : num  5087 463 5802 4667 3843 ...
 $ P0080010  : num  4703 445 5296 4347 3577 ...
 $ P0080011  : num  2181 181 2242 2165 1817 ...
 $ P0080012  : num  843 160 467 1094 588 ...
 $ P0080013  : num  797 36 847 397 576 ...
 $ P0080014  : num  55 2 16 45 27 36 134 19 8 13 ...
 $ P0080015  : num  318 47 348 230 226 ...
 $ P0080016  : num  147 5 198 202 91 175 523 31 12 19 ...
 $ P0080017  : num  68 0 111 42 44 65 348 10 10 9 ...
 $ P0080018  : num  29 0 64 16 8 13 53 2 4 5 ...
 $ P0080019  : num  120 3 271 77 93 160 953 21 7 15 ...
 $ P0080020  : num  18 1 72 22 9 39 199 15 2 3 ...
 $ P0080021  : num  2 0 3 2 3 8 12 4 1 1 ...
 $ P0080022  : num  8 2 107 14 13 41 151 3 1 2 ...
 $ P0080023  : num  26 1 28 10 11 15 70 8 4 9 ...
 $ P0080024  : num  80 7 420 29 62 51 620 9 4 7 ...
 $ P0080025  : num  11 0 102 2 9 18 95 5 1 2 ...
 $ P0080026  : num  354 18 489 298 238 360 819 54 46 82 ...
 $ P0080027  : num  202 7 236 191 135 225 380 33 24 50 ...
 $ P0080028  : num  67 1 36 30 17 33 66 1 4 7 ...
 $ P0080029  : num  5 2 3 2 7 3 10 2 4 0 ...
 $ P0080030  : num  23 5 56 24 29 36 87 3 2 3 ...
 $ P0080031  : num  12 0 7 16 18 7 29 2 6 4 ...
 $ P0080032  : num  0 0 4 1 2 2 7 0 1 0 ...
 $ P0080033  : num  5 1 3 5 4 11 20 4 1 2 ...
 $ P0080034  : num  16 0 11 12 3 17 20 3 1 7 ...
 $ P0080035  : num  9 1 7 2 4 8 91 1 2 4 ...
 $ P0080036  : num  2 0 5 2 2 2 21 1 1 0 ...
 $ P0080037  : num  2 0 8 5 1 7 20 0 0 3 ...
  [list output truncated]

df[df$NAME10 == 'Broome',]['fips']

Question 1(e)

Compute descriptive statistics of the population column (P0010001), including total, minimum, maximum, mean, median, and skewness.

if(!require(moments)) install.packages('moments')
library(moments)
population <- df$P0010001
total <- sum(population)
minimum <- min(population)
maximun <- max(population)
mean <- mean(population)
median <- median(population)
skewness <- skewness(population)
total

[1] 19378102

minimum

[1] 4836

maximun

[1] 2504700

mean

[1] 312550

median

[1] 91301

skewness

[1] 2.579801

Question 1(f)

Create a histogram and a boxplot of the population.

hist(population)

boxplot(population)

Question 2

Use R package UScensus2010tract to complete the following tasks: (20 pt.)

if(!require('UScensus2010tract')) install.tract('osx')
library(UScensus2010tract)

Question 2(a)

Plot a map of New York census tracts using the plot function.

data(new_york.tract10)
shp <- new_york.tract10
plot(shp)

Question 2(b)

Compute the total population based on census tracts.

df <- shp@data
sum(df$P0010001)

[1] 19378102

df

Question 2(c)

Select all census tracts in Broome County and plot the map.

plot(shp[df$county == '007',])

Question 2(d)

What’s the total population of Broome County?

broome <- shp[df$county == '007',]
broomedf <- broome@data
sum(broomedf$P0010001)

[1] 200600

Question 2(e)

Create a histogram and a boxplot of population based on census tracts of Broome County.

hist(broomedf$P0010001)

boxplot(broomedf$P0010001)

Question 2(f)

Select the first five columns of the shapefile of Broome County census tract; add a new population ratio column (= census tract population / county population); save the new shapefile to the hard drive.

broome_5 <- broome[,1:5]
df <- broome_5@data
df$popratio <- df$P0010001 / sum(broomedf$P0010001)
broome_5@data <- df
if(!require(rgdal)) install.packages('rgdal')
library(rgdal)
getwd()

[1] "C:/Users/hwang164/Documents/GEOG-533/Labs/Lab_09"

writeOGR(broome_5,dsn = '.', layer = "broome_5", driver = 'ESRI Shapefile')

Error in writeOGR(broome_5, dsn = ".", layer = "broome_5", driver = "ESRI Shapefile") : 
  layer exists, use a new layer name

Question 3

Use R packages raster and ncdf4 to complete the following tasks: (20 pt.)

if(!require('raster')) install.packages('raster')
library(raster)
if(!require('ncdf4')) install.packages('ncdf4')

Loading required package: ncdf4
there is no package called 㤼㸱ncdf4㤼㸲Installing package into 㤼㸱C:/Users/hwang164/Documents/R/win-library/3.5㤼㸲
(as 㤼㸱lib㤼㸲 is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/3.5/ncdf4_1.16.1.zip'
Content type 'application/zip' length 2598248 bytes (2.5 MB)
downloaded 2.5 MB

package ‘ncdf4’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\hwang164\AppData\Local\Temp\RtmpOegO0P\downloaded_packages

library(ncdf4)

package 㤼㸱ncdf4㤼㸲 was built under R version 3.5.3

Question 3(a)

Load the multi-band raster dataset “NDVI.nc” into R.

ndvi.rb <- brick('NDVI.nc')

Question 3(b)

Get the basic information about the dataset, including the number of rows, columns, cells, and bands; spatial resolution, extent, bounding box, and projection.

nrow(ndvi.rb)

[1] 360

ncol(ndvi.rb)

[1] 720

ncell(ndvi.rb)

[1] 259200

nlayers(ndvi.rb)

[1] 36

res(ndvi.rb)

[1] 0.5 0.5

extent(ndvi.rb)

class       : Extent 
xmin        : -180 
xmax        : 180 
ymin        : -90 
ymax        : 90 

bbox(ndvi.rb)

    min max
s1 -180 180
s2  -90  90

projection(ndvi.rb)

[1] "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

Question 3(c)

Aggregate all bands to generate a mean NDVI raster and a maximum NDVI raster; save the two new raster datasets to the hard drive.

mean <- mean(ndvi.rb)
max <- max(ndvi.rb)
plot(mean)

plot(max)

writeRaster(mean, filename = 'mean_ndvi')

class       : RasterLayer 
dimensions  : 360, 720, 259200  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : C:/Users/hwang164/Documents/GEOG-533/Labs/Lab_09/mean_ndvi.grd 
names       : layer 
values      : 0, 0.8537528  (min, max)

writeRaster(max, filename = 'max_ndvi')

class       : RasterLayer 
dimensions  : 360, 720, 259200  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : C:/Users/hwang164/Documents/GEOG-533/Labs/Lab_09/max_ndvi.grd 
names       : layer 
values      : 0, 0.9681  (min, max)

Question 3(d)

Plot the maps of monthly NDVI of the year 2001

Question 3(e)

Create histograms of monthly NDVI of the year 2001.

hist(ndvi2001)

size changed to n because it cannot be larger than n when replace is FALSE39% of the raster cells were used (of which 42% were NA). 58093 values used.

Question 3(f)

Plot the NDVI map of July 2000; click any location with data on the map and retrieve the NDVI time series for all years; plot the NDVI time series of the selected location.

values <- click(ndvi.rb, n = 1, xy = FALSE)

Error in graphics::locator(1, type, ...) : 
  plot.new has not been called yet